subroutine agp_local_energy_adjoints_aa_aaaa_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
       RT, & ! (in)    unoccupied row matrix times inverse of occ pairing matrix
      RTt, & ! (in)    transpose of RT
      Tia, & ! (in)    occ-unocc slice of the one-electron integrals
    Viajb, & ! (in)    occ-unocc-occ-unocc slice of the two-electron integrals
        E, & ! (out)   energy contribution
      RTA)   ! (out)   adjoint for RT

implicit none

integer,      intent(in)    ::    ni
integer,      intent(in)    ::    na
real(kind=8), intent(in)    ::    RT(na,ni)
real(kind=8), intent(in)    ::    RTt(ni,na)
real(kind=8), intent(in)    ::    Tia(ni,na)
real(kind=8), intent(in)    ::    Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::    E
real(kind=8), intent(out)   ::    RTA(ni,na)

integer :: i, a, j, b

real(kind=8) :: tempX, tempC, tempE, tempJ

! initialize outputs to zero
E = 0.00d+00
RTA = 0.00d+00

! one electron integrals
do a = 1,na
do i = 1,ni

  tempJ = Tia(i,a)

  RTA(i,a) = RTA(i,a) + tempJ

  E = E + RTt(i,a) * tempJ

enddo
enddo

! loop over two electron integrals
do b = 1,na
  do j = 1,ni
    do a = 1,na
      do i = 1,ni

        ! get tei element
        tempJ = Viajb(i,a,j,b)

        ! compute coulomb and exchange products with RT element
        tempC = tempJ * RTt(i,a)
        tempX = tempJ * RT(a,j)

        ! add RT adjoint coulomb and exchange contributions
        RTA(j,b) = RTA(j,b) + tempC
        RTA(i,b) = RTA(i,b) - tempX

        ! compute energy contribution
        tempE = tempC * RTt(j,b) - tempX * RTt(i,b)

        ! add energy contribution
        E = E + 0.50d+00 * tempE

      enddo
    enddo
  enddo
enddo

end subroutine agp_local_energy_adjoints_aa_aaaa_f




subroutine agp_local_energy_adjoints_aabb_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
      RTt, & ! (in)    transpose of contraction of unoccupied rows with inverse matrix
       TC, & ! (in)    contraction of unoccupied columns with inverse matrix
      RTC, & ! (in)    contraction of RT with unoccupied columns
     IPMt, & ! (in)    transpose of the inverse pairing matrix
      UPM, & ! (in)    the unoccupied pairing matrix
    Viajb, & ! (in)    two electron integrals in 1122 order
        E, & ! (out)   energy contribution
      VRT, & ! (out)   contraction of integrals with the RT matrix
      VTC, & ! (out)   contraction of integrals with the TC matrix
       VT, & ! (out)   contraction of integrals with inverse pairing matrix
    VARTC)   ! (out)   contraction of integrals with ( UPM - RTC )

implicit none

integer,      intent(in)    ::   ni
integer,      intent(in)    ::   na
real(kind=8), intent(in)    ::   RTt(ni,na)
real(kind=8), intent(in)    ::   TC(ni,na)
real(kind=8), intent(in)    ::   RTC(na,na)
real(kind=8), intent(in)    ::   IPMt(ni,ni)
real(kind=8), intent(in)    ::   UPM(na,na)
real(kind=8), intent(in)    ::   Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::   E
real(kind=8), intent(out)   ::   VRT(ni,na)
real(kind=8), intent(out)   ::   VTC(ni,na)
real(kind=8), intent(out)   ::   VT(na,na)
real(kind=8), intent(out)   ::   VARTC(ni,ni)

integer :: i, a, j, b

real(kind=8) :: tempE, tempV, tempARTC, tempQ

E = 0.00d+00
VRT = 0.00d+00
VTC = 0.00d+00
VT = 0.00d+00
VARTC = 0.00d+00

do b = 1,na

  do j = 1,ni

    do a = 1,na

      tempARTC = UPM(a,b) - RTC(a,b)

      do i = 1,ni

        ! get integral
        tempV = Viajb(i,a,j,b)

        ! process terms involving integral times RT
        tempQ = tempV * RTt(i,a)
        tempE = tempQ * TC(j,b)
        VRT(j,b) = VRT(j,b) + tempQ

        ! process terms involving integral times TC
        VTC(i,a) = VTC(i,a) + tempV * TC(j,b)

        ! process terms involving integral times inverse pairing matrix  (a.k.a. T)
        tempQ = tempV * IPMt(i,j)
        tempE = tempE + tempQ * tempARTC
        VT(a,b) = VT(a,b) + tempQ

        ! process terms involving integral times ARTC term
        VARTC(i,j) = VARTC(i,j) + tempV * tempARTC

        ! add energy contribution
        E = E + tempE

      enddo
    enddo
  enddo
enddo

end subroutine agp_local_energy_adjoints_aabb_f
